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Abstract 

A new method for the computation of flows at the micrometer scale is presented. 
It is based on the recently introduced minimal entropic kinetic models. Both the 
thermal and isothermal families of minimal models are presented, and the simplest 
isothermal entropic lattice Bhatnagar-Gross-Krook (ELBGK) is studied in detail 
in order to quantify its relevance for microflow simulations. ELBGK is equipped 
with boundary conditions which are derived from molecular models (diffusive wall). 
A map of three-dimensional kinetic equations onto two-dimensional models is es- 
tablished which enables two-dimensional simulations of quasi- two-dimensional flows. 
The ELBGK model is studied extensively in the simulation of the two-dimensional 
Poiseuille channel flow. Results are compared to known analytical and numerical 
studies of this flow in the setting of the Bhatnagar-Gross-Krook model. The ELBGK 
is in quantitative agreement with analytical results in the domain of weak rarefaction 
(characterized by Knudsen number Kn, the ratio of mean free path to the hydrody- 
namic scale), up to Kn ~ 0.01, which is the domain of many practical microflows. 
Moreover, the results qualitatively agree throughout the entire Knudsen number 
range, demonstrating Knudsen's minimum for the mass flow rate at moderate values 
of Kn, as well as the logarithmic scaling at large Kn. The present results indicate 
that ELBM can complement or even replace computationally expensive microscopic 
simulation techniques such as kinetic Monte Carlo and/or molecular dynamics for 
low Mach and low Knudsen number hydrodynamics pertinent to microflows. 

1 Introduction 

Gas flows at the micrometer scale constitute a major portion of contemporary fluid dy- 
namics of engineering interest. Because of its relevance to the engineering of micro electro- 
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mechanical systems (MEMS), the branch of computational fluid dynamics focused on micro 
scale phenomena is often called "microfluidics" jU |2] • 

Microflows are characterized by the Knudsen number, Kn, which is deflned as the 
ratio of the mean free path of molecules A and the characteristic scale L of variation of 
hydrodynamic flelds (density, momentum, and energy). For typical flows in microdevices, 
Kn ~ \/L varies from Kn ^ 1 (almost- continuum flows) to Kn ~ 1 (weakly rarefled 
flows). Another characteristic property of microflows is that they are highly subsonic, that 
is, the characteristic flow speed is much smaller than the speed of sound. This feature is 
characterized by the Mach number. Ma ~ m/cs, where u is the characteristic flow speed, 
and Cs is the (isentropic) speed of sound. Thus, for microflows. Ma <^ 1. To be more 
speciflc, typical flow velocities are about 0.2 m/s, corresponding to Ma ~ 10~^, while 
values of the Knudsen number range between 10~^ < Kn < 10^^. Finally, in the majority 
of applications, microflows are quasi-two-dimensional. 

Theoretical studies of gas flows at flnite Knudsen number have begun several decades 
ago in the realm of the Boltzmann kinetic equation. To that end, we mention pioneering 
contributions by Cercignani, Sone, and others Olll- These studies focused on obtaining 
either exact solutions of the stationary Boltzmann kinetic equation, or other model kinetic 
equations in relatively simple geometries (most often, inflnite or semi-inflnite rectangular 
ducts), or asymptotic expansions of these solutions. 

While analytical solutions are important for a qualitative understanding of microflows, 
and also for the validation of numerical schemes, they certainly do not cover all the needs 
of computational fluid dynamics of practical interest. At present, two CFD strategies for 
microflows are well established. 

• Equations of continuous fluid mechanics with slip boundary conditions. The simplest 
semi-phenomenological observation about microflows is the break down of the no-slip 
boundary condition of fluid mechanics with increasing Knudsen number. Since mi- 
croflows are highly subsonic, this leads to the simplest family of models, equations of 
incompressible or compressible fluid dynamics supplemented by slip velocity bound- 
ary conditions (a review can be found in [2]). This approach, although widely used 
at the early days of microfluidics, remains phenomenological. Moreover, it fails to 
predict phenomena such as non-trivial pressure and temperature proflles observed by 
more microscopic approaches. 

• Direct simulation of the Boltzmann kinetic equation. On the other extreme, it is 
possible to resort to a fully microscopic picture of collisions, and to use a molecular 
dynamics approach or a simplifled version thereof - the Direct Simulation Monte 
Carlo method of Bird (DSMC) DSMC is sometimes heralded as the method 
of choice for simulation of the Boltzmann equation, and it has indeed proven to be 
robust in supersonic, highly compressible flows with strong shock waves. However, 
the highly subsonic flows at small to moderate Knudsen number is not a "natural" 
domain for the DSMC simulations where it becomes computationally intensive jB]. 

Since semi-phenomenological computations are not reliable, and the fully microscopic 
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treatment is not feasible, the approach to CFD of microflows must rely on reduced models 
of the Boltzmann equation. Two classical routes of reducing the kinetic equations are 
well known, the Chapman-Enskog method and Grad's moment method (for a modern 
summary and extensions of these methods see, for example, [Zj). The Chapman-Enskog 
method extends the hydrodynamic description (compressible Navier-Stokes equations) to 
finite Kn in the form of a Taylor series, leading to hydrodynamic equations of increasingly 
higher order in the spatial derivatives (Burnett's hydrodynamics). Grad's method extends 
the hydrodynamic equations to a closed set of equations including higher-order moments 
(fluxes) as independent variables. Both methods are well suited for theoretical studies 
of microflows. In particular, as was already noted by Grad jHj, moment equations are 
especially well suited for low Mach number flows. 

However, applications of Grad's moment equations or of Burnett's hydrodynamics (or 
of existing extensions and generalizations thereof) to CFD of microflows are limited at 
present because of several reasons. The most severe difficulty is in formulating the boundary 
conditions at the reduced level. Although some studies of boundary conditions for moment 
systems were initiated recently , this problem is far from solved. The crucial importance 
of the boundary condition for microflows is actually expected. Indeed, as the rarefaction 
is increasing with Kn, the contribution of the bulk collisions becomes less significant as 
compared to the collisions with the boundaries, and thus the realistic modelling of the 
boundary conditions becomes increasingly important. 

In this paper we set up a novel approach to the CFD of microfiows. It is based on the 
recently developed Entropic Lattice Boltzmann Method (ELBM) [101 CH 1121 113 El 113 CHI 
El El CHI- The choice of the ELBM for microfiows is motivated by two reasons: 

• ELBM is an unconditionally stable simulation method for fiows at low Mach numbers. 

• In contrast to Grad's method, ELBM is much more compliant with the boundary con- 
ditions. Recently, an appropriate boundary condition for the ELBM was found upon 
a discretization of the diffusive wall boundary condition of the Boltzmann equation 
jl5j . This boundary condition was also rediscovered in [20] , where ELBM simulations 
were tested against molecular dynamic simulations with a good agreement. 

It should be mentioned that the predecessor of ELBM, the lattice Boltzmann method 
(LBM) [21], was employed several times for microflow simulations [221 123 El I2S]- The 
motivation of most of these works was the velocity slip observed in the LBM simulations 
using the so-called bounce-back boundary condition [221 123 1211 12S] ■ Hovewer, since the 
bounce-back boundary condition is completely artificial, the results are questionable fK)\ . 

The outline of the present paper is as follows: In section |21 for the sake of completeness, 
we present the general description of isothermal and thermal entropic lattice Boltzmann 
models. Then, we describe the ELBM setup for the simplest situation of isothermal fiows, 
the entropic lattice Bhatnagar-Gross-Krook model (ELBGK). The case of isothermal mod- 
els is important in itself since it allows to study nontrivial fiow phenomena such as velocity 
slip at the wall, and the well-known Knudsen minimum of the mass fiowrate in pressure- 
driven channel fiows. Prediction of the Knudsen minimum is a classical benchmark problem 
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for the simulation of microflows. In sectional we describe in detail the implementation of 
the diffusive boundary conditions for the ELBM for two-dimensional simulations. A sep- 
arate section E] is devoted to the question of how to simulate quasi-two-dimensional flows 
with two-dimensional models, and how to map the parameters of the model onto experi- 
mental data and more microscopic simulations. In sectional we present simulation results 
for the quasi-two-dimensional Poiseuille flow, and discuss Knudsen's minimum, comparing 
results with known asymptotic and analytic solutions to the Boltzmann kinetic equation. 
Results are discussed in section IHl where we define the domain of validity of ELBM for 
microflows. We also suggest a straightforward application of ELBM results to accelerate 
more microscopic simulation approaches like the DSMC method. 

2 Minimal kinetic models 

We start with a generic discrete velocity kinetic model. Let fi{x,t) be populations of 
the D- dimensional discrete velocities Cj, i = 1, ...,nd, at position x and time t. The 
hydrodynamic fields are the linear functions of the populations, namely 

"d 

d, c^,}fi = {p, pu, pDT + pn^}, (1) 

2 = 1 

where p is the mass density of the fluid, pu is the dimensional momentum density vec- 
tor, and e = pDT + pv? is the energy density. In the case of isothermal simulations, the 
set of independent hydrodynamic fields contains only the mass and momentum densities. 
It is convenient to introduce rid-dimensional population vectors f , and the standard scalar 
product, (f|g) = T^ti^iUi- Foi' example, for almost-incompressible hydrodynamics (leav- 
ing out the energy conservation), the locally conserved density and momentum density 
fields are written as 

(l|f)=p, (c,|f)=pn„. (2) 

Here 1 = {1}"^]^, Cq, = {cja}"ii, a = 1, . . . , where D is the spatial dimension. 

The construction of the kinetic simulation scheme begins with finding a convex function 
of populations H (entropy function), which satisfies the following condition: If f°''(p, u) 
(local equilibrium) minimizes H subject to the hydrodynamic constraints (equations (0) or 
(121)), then f''^ also satisfies certain restrictions on the higher-order moments. For example, 
the equilibrium stress tensor must respect the Galilean invariance, 

CiaCipfi\p, U) = pclSa/3 + pUaUp. (3) 

1=1 

The corresponding entropy functions for the isothermal and the thermal models were found 
in fl^ [W\ and are given below (see section l!^.U.4l and Table For the time being, 
assume that the convex function H is given. 
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The next step is to obtain the set of kinetic equations, 

dtfi + Ciadafi = Ai. (4) 

Let mi, ... , be the na-dimensional vectors of locally conserved fields. Mi — (mi|f), 
i = 1, . . . ,nc, He < Ud- The nd-dimensional vector function A (collision integral), must 
satisfy the conditions: 

(mj| A) = 0, i = 1, . . . , ric (local conservation laws), 

(V-?J|A) < (entropy production inequality), 

where V/J is the row-vector of partial derivatives dH/dfi. Moreover, the local equilibrium 
vector f'^'^ must be the only zero point of A, that is, A(f'''^) = 0, and, finally, P'^ must be 
the only zero point of the local entropy production, a(f°'^) — 0. CoUision integrals which 
satisfies all these requirements are called admissible. Let us discuss several possibilities of 
constructing admissible collision integrals. 

2.0.1 BGK model 

Suppose that the entropy function H is known. If, in addition, the local equilibrium 
is also known as an explicit function of the locally conserved variables (or some reliable 
approximation of this function is known), the simplest option is to use the Bhatnagar- 
Gross-Krook (BGK) model. In the case of isothermal hydrodynamics, for example, we 
write 

A = -i(f-r^(p(f),n(f))). (5) 

The BGK collision operator is sufficient for many applications. However, it becomes advan- 
tageous only if the local equilibrium is known in a closed form. Unfortunately, often only 
the entropy function is known but not its minimizer. For these cases one should construct 
collision integrals based solely on the knowledge of the entropy function. We present here 
two particular realizations of the collision integral based on the knowledge of the entropy 
function only. 

2.0.2 Quasi-chemical model 

For a generic case of ric locally conserved fields, let g^, s = 1, . . . , na — be a basis of the 
subspace orthogonal (in the standard scalar product) to the vectors of the conservation 
laws. For each vector g^, we define a decomposition g^ = g+ — gj, where all components 
of vectors g^ are nonnegative, and if gf^ ^ 0, then gf^ — 0. Let us consider the collision 
integral of the form: 

^ = E ^^^s {exp ((Vi/ |g;)) - exp ((Vi/ |g+)) } . (6) 
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Here > 0. By construction, the collision integral © is admissible. If the entropy 
function is Boltzmann-like, and the components of the vectors are integers, the collision 
integral assumes the familiar Boltzmann-like form. An example of such a collision term 
for the D2(59-discrete velocity model is described in Ref. j26j . 

2.0.3 Single relaxation time gradient model 

The BGK collision integral (0) has the following important property: the linearization of 
the operator (0) around the local equilibrium point has a very simple spectrum {0, — 1/r}, 
where is the ric-times degenerate eigenvalue corresponding to the conservation laws, while 
the non-zero eigenvalue corresponds to the rest of the (kinetic) eigenvectors. Nonlinear 
collision operators which have this property of their linearizations at equilibrium are called 
single relaxation time models (SRTM). They play an important role in modelling because 
they allow for the simplest identification of transport coefficients. 

The SRTM, based on the given entropy function if, is constructed as follows (single 
relaxation time gradient model, SRTGM). For the system with Uc local conservation laws, 
let e^, s = 1, . . . , nd — ric, be an orthonormal basis in the kinetic subspace, (mj|es) = 0, 
and (e^lep) = 5sp. Then the single relaxation time gradient model is 

A = — 5^ e,K,,{i){e,\VH), (7) 

s,p=l 

where Ksp are elements of a positive definite {rid — n^) x (nd — Uc) matrix 

K{f) = C~\f), (8) 
C,p(f) = (e,|VVii(f)|e,). 

Here, W/f(f) is the x matrix of second derivatives, d'^H/dfidfj. Linearization of 
the collision integral at equilibrium has the form, 

L = — V e,e„ (9) 

which is obviously single relaxation time. Use of the SRTGM instead of the BGK model 
results in the same hydrodynamics even when the local equilibrium is not known in a closed 
form. Further details of this model can be found in Ref. jl4j . 

2.0.4 if-functions of minimal kinetic models 

The Boltzmann entropy function written in terms of the one-particle distribution function 
f{x, c) is H = J / In / dc, where c is the continuous velocity. Close to the global (reference) 
equilibrium, this integral can be approximated by using the Gauss-Hermite quadrature 
with the weight 

W^ = (2 7rro)(^/2)exp(-cV(2ro)). 
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Table 1: Minimal kinetic models. Column 1: Order of Hermite velocity polynomial used 
to evaluate the Gauss-Hermite quadrature; Column 2: Locally conserved (hydrodynamic) 
fields; Column 3: Discrete velocities for D = 1 (zeroes of the corresponding Hermite 
polynomials). For D > 1, discrete velocities are all possible tensor products of the one- 
dimensional velocities in each component direction; Column 4: Weights in the entropy 
formula ()10|) . corresponding to the discrete velocities of the Column 3. For D > 1, the 
weights of the discrete velocities are products of corresponding one-dimensional weights; 
Column 5: Macroscopic equations for the fields of Column 2 recovered in the hydrodynamic 
limit of the model. 



1. Order 


2. Fields 


3. Velocities 


4. Weights 


5. Hydrodynamic limit 
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1/2 
1/2 
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p, pu 





2/3 
1/6 
1/6 


Isothermal Navier-Stokes 


4 


p, pu, e 


V^ + VqVTo 
-V^ + VgVTo 


l/[4(3 - v^)] 
l/[4(3 - v^)] 
l/[4(3 + VQ)] 
l/[4(3 + v^)] 


Thermal Navier-Stokes 



Here D is the spatial dimension, Tq is the reference temperature, while the particles mass 
and Boltzmann's constant ks are set equal to one. This gives the entropy functions of the 
discrete- velocity models [2^11131111, 

H^tf'^'iQ- (10) 

i=l ^ ' 

Here, Wi is the weight associated with the i-th discrete velocity c, (zeroes of the Hermite 
polynomials). The discrete- velocity distribution functions (populations) fi{x) are related 
to the values of the continuous distribution function at the nodes of the quadrature by the 
formula, 

/,(a.) = m,(27rTo)(^/2)exp(c^/(2To))/(x,Q). 

The entropy functions ()10|) for various {u'j, c, } are the only input needed for the construc- 
tion of minimal kinetic models. 

With the increase of the order of the Hermite polynomials used in evaluation of the 
quadrature (fTUI) . a better approximation to the hydrodynamics is obtained. The first few 
models of this sequence are represented in Tabled 
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2.1 Entropic lattice Boltzmann method 



If the set of discrete velocities forms the hnks of a Bravais lattice (with possibly several 
sub- lattices) , then the discretization of the discrete velocity kinetic equations in time and 
space is particularly simple, and leads to the entropic lattice Boltzmann scheme. This 
happens in the important case of the isothermal hydrodynamics. The equation of the 
entropic lattice Boltzmann scheme reads 

fi{x + c,6t, t + 6t)- fi{x, t) = (3a{i{x, t))A,{i{x, t)), (11) 

where 6t is the discretization time step, and (3 G [0, 1] is a fixed parameter which matches 
the viscosity coefficient in the long-time large-scale dynamics of the kinetic scheme flllj) . 
The function a of the population vector defines the maximal over-relaxation of the scheme, 
and is found from the entropy condition, 

H{f{x, t) + aA(f (a;, t)) = H{f{x, t)). (12) 

The nontrivial root of this equation is found for populations at each lattice site. Equation 
()12|) ensures the discrete-time i7-theorem, and is required in order to stabilize the scheme 
if the relaxation parameter (3 is close to one. We note in passing that the latter limit 
is of particular importance in the applications of the entropic lattice Boltzmann method 
to hydrodynamics because it corresponds to vanishing viscosity, and hence to numerically 
stable simulations of very high Reynolds number flows. 



2.2 Entropic lattice BGK method (ELBGK) 

An important simpliflcation occurs in the case of the isothermal simulations when the 
entropy function is constructed using third-order Hermite polynomials (see Table CJ: the 
local equilibrium population vector can be obtained in closed form JHl- This enables the 
simplest entropic scheme - the entropic lattice BGK model - for simulations of isothermal 
hydrodynamics. We present this model in dimensionless lattice units. 

Let D be the spatial dimension. For D = 1, the three discrete velocities are 

c = {-1,0,1}. (13) 

For D > 1, the discrete velocities are tensor products of the discrete velocities of these one- 
dimensional velocities. Thus, we have the 9-velocity model for D = 2 and the 27-velocity 
model for D = 3. The H function is Boltzmann-like: 

i=l ^ ' 

The weights Wi are associated with the corresponding discrete velocity Cj. For D = 1, the 
three-dimensional vector of the weights corresponding to the velocities ()13|) is 
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For D > 1, the weights are constructed by multiplying the weights associated with each 
component direction. 

The local equilibrium minimizes the if-function (11 Oj) subject to the fixed density and 
momentum, 

3^ 3^ 

^fi = XI ^ a = 1, • • • 5 ^- (16) 

i=l i=l 

The explicit solution to this minimization problem reads, 



Note that the exponent, Cia, in (fTTj) takes the values ±1, and only, and the speed of 
sound, Cs, in this model is equal to 1/a/3. The factorization of the local equilibrium (fT7|) 
over spatial components is quite remarkable, and resembles the familiar property of the 
local Maxwellians. 

The entropic lattice BGK model for the local equilibrium (fTTj) reads, 

fi{x + c,6t, t + 6t)- U{x, t) = -Pa{f,{x, t) - /r(p(f (a;, t)), n(f (a;, t))). (18) 

The parameter j3 is related to the relaxation time r of the BGK model (0) by the formula, 

and the value of the over-relaxation parameter a is computed at each lattice site from the 
entropy estimate, 

iy(f-a(f-f"^(f))) = i7(f). (20) 

In the hydrodynamic limit, the model (jl8|) reconstructs the Navier-Stokes equations with 
the viscosity 

^i = pclr = pcl5t - • (21) 
The zero- viscosity limit corresponds to /? — 1. 

3 Wall boundary conditions 

The boundary (a solid wall) OR is specified at any point x G OR by the inward unit normal 
e, the wall temperature T^, and the wall velocity u^. The simplest boundary condition 
for the minimal kinetic models is obtained upon evaluation of the diffusive wall boundary 
condition for the Boltzmann equation [3^ with the help of the Gauss-Hermite quadrature 
[T^ Uni 120] • The essence of the diffusive boundary condition is that particles loose their 
memory of the incoming direction after reaching the wall. Once a particle reaches the wall, 
it gets redistributed in a way consistent with the mass-balance and normal-flux conditions. 
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Further, the boundary condition must also satisfy the condition of detailed balance: if the 
incoming populations are at equilibrium (corresponding to the wall- velocity), the outgoing 
populations are also at equilibrium (corresponding to the wall- velocity). 

For the purpose of simulations below, let us consider the case when the wall normal, 
e, (pointing towards the fluid) is in the positive y direction. The lattice for the for the 
9- velocity isothermal model is depicted in Fig. ^ 
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Figure 1: Schematic diagram for the situation near a flat wall, when the wall normal, e, 
(pointing towards the fluid) is in the positive y direction. 

For this particular case, the boundary update rules for incoming and grazing popula- 
tions on a two-dimensional lattice are: 

fo{x,y,t + 6t) = f*{x,y,t), 
fi{x,y,t + 6t) = /o(x - c6t,y,t), 
fsix, y, t + 6t) = /o (x + c6t, y, t), 

U{x, y,t + 6t) = ^ [f:{x, y + c6t, t) + /;(x, y, t)] , (22) 
/7(x, y,t + 5t) = ]^ [f;{x + c5t, y + c5t, t) + /^(x, y, t)] , 
/8(x, t + 5t) = i [f*{x - c6t, y + c6t, t) + /§ (x, y, t)] , 
where /* denotes post-collision populations, and the update rules for outgoing populations 



are: 



f-^( ./4(x,t + 5t) + /7(x,t + 5t) + /8(x,t-f5t) 



/6(x,t + 5t) = /|'l(p,U 



wall J 



fTip, Uwall) + JTiP^ U^all) + ftiP^ "wall) 
ft^,,,,^ .eq. , /4(X, t + 6t) + /7(X, t + 6t) + fs{^, t + 5t) 

/2 (P, U„all) + /s (P, Uwall) + /e (P, Uwall) 
^ /4(X, t + 5t) + /7(X, t + 5t) + /8(X, t + 5t) 

tTi.P^ U^all) + tT{P^ Uwall) + fl'^ip, U^all) ' 
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4 Simulation of quasi-two-dimensional flows 
with two-dimensional kinetic models 



As mentioned in the introduction, many microflows of engineering interest can be consid- 
ered as quasi-two-dimensionaL This means that averaged quantities such as flow velocity 
and density do not depend appreciably on the third spatial direction. Thus, it is tempting 
to use two-dimensional kinetic models in simulations of such flows. However, care should 
be taken in order to map correctly the results of two-dimensional simulations onto experi- 
mental data or molecular dynamics simulations. Indeed, molecular motion remains three- 
dimensional in spite of the fact that some averages can be considered two-dimensional. In 
the DSMC simulations of two-dimensional flows, for example, collisions of the particles 
are always treated as three-dimensional. The two-dimensional kinetic models therefore 
must be considered as a computational device which uses fictitious particles moving in 
two dimensions in order to mimic quasi-two-dimensional flows of particles moving in three- 
dimensions. 

The mapping of the parameters of the three-dimensional kinetic equation onto the 
two-dimensional lattice Boltzmann scheme is done in two steps: 

• Map the continuous 3D kinetic equation onto the 3D discrete velocity model. 

• Map the 3D discrete velocity model onto the 2D velocity model. 

In the case considered in this paper, the 3D continuous kinetic model is the BGK model 
|3j, which contains the relaxation parameter tbgk, 

dtf + cM = —if-f''), (24) 

Tbgk 

where is the local Maxwell distribution function. In this subsection we shall explicitly 
indicate all the functions and parameters related to the continuous BGK model with the 
subscript in order to distinguish them from the lattice counterparts. 

The viscosity coefficient of the BGK model ()24j) is related to the relaxation time tbgk 
as follows: 

fJ'BGK = PbgkTbgkTbgk- (25) 

In the isothermal discrete velocity BGK model of section |2l the viscosity is expressed 
not as a function of temperature but of sound speed Cg, 

fi = pc^T, (26) 

where Cg = 1/v^ in lattice units used in simulations. Therefore, in the first step of the 
above procedure, we map the parameters of the continuous 3D BGK equation onto the 3D 
discrete velocity model using the relation for the speed of sound, Cg bgk = Vt^bck, where 
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7 is the adiabatic exponent. For an ideal gas, 7 = 5/3. Thus, the first step is accomphshed 
by the relation. 



Cs BGK — Cs, 




This formula establishes the relation between the relaxation parameter of the 3D contin- 
uous BGK equation (j24|) and the 3D isothermal discrete velocity model of section 
Note that the sound speed of the thermal model is made equal to the sound speed of the 
isothermal model by the relation ()27p. 

At the second step, we map the three-dimensional 27-velocity isothermal model onto the 
two-dimensional 9-velocity model. Note that the sound speed in both models is identical 
(and equals to V/cb^o in dimensional units). The mapping is done by populating at the 
equilibrium the links of the 27-velocity lattice in the direction orthogonal to the fixed 
plane containing the links of the 9-velocity sub-lattice. This amounts to the following 
recomputation of the three-dimensional density p^n in terms of the two-dimensional density 

P2D- 

3 

P3D = ^P2D- (28) 

This formula enables the computation of the effective three-dimensional density in 
terms of the two-dimensional density used in the simulations with the 9-velocity model. 
Specification of density is an important part of the simulation since it is used to define 
data as prescribed pressure drop at the inlet and outlet of pipes etc. 

The formulas collected in this section, together with the viscosity formula of the 
fully discretized entropic lattice Boltzmann method ()2H1 . enable the comparison of two- 
dimensional simulation results with the results obtained by microscopic simulations with 
different collision models (the hard sphere model, for example). The choice of the model is 
needed to identify the mean free path which is is model-dependent, especially at moderate 
values of Knudsen number. 

5 Plane Poiseuille flow 

Plane Poiseuille flow is one of the most studied benchmarks on gas dynamics. The gas 
moves between two parallel plates driven by a flxed pressure difference between the inlet 
and outlet. It is well known that for this setup the flow rate through a cross-section of the 
pipe exhibits a minimum [2ZIIS]- In fact, one of the major achievements in the early days 
of kinetic theory was the prediction of a minimum of the mass flow rate as a function of 
the Knudsen number for Kn ~ 1. 

We simulate the two-dimensional flow in a rectangular duct of length L along the 
streamwise direction (x) and width H <^ L along the wall-normal {y) direction. The flow 
is driven by a flxed pressure difference AP = Pm — -Pout) where Pin and Pout are the pressure 
at the inlet and outlet of the duct, respectively. In the subsequent analysis, we shall follow 
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the convention used by Cercignani, where the Knudsen number for the BGK model is 
defined as (continuous BGK units): 



Kn = (29) 

where the pressure Pq is defined as the mean of the inlet and outlet pressures, Pq = 
(Pin + Pout)/2. In the hydrodynamic limit, this results in the well known Hagen-Poiseuille 
parabolic velocity profile: 

uiy) = Uo(^\- I') , (30) 

where the amplitude of the flow for a two-dimensional duct is Uq = H'^AP/{2fiL). From 
the analysis of the Boltzmann-BGK equation it is known that the dimensionless flow rate 

Q = -—— / u{y)dy, (31) 

has a low-Knudsen (Kn ^ 1) asymptotic: 

6Kn 

with s = 1.015 (see 0), which is equivalent to solving the Navier-Stokes equation with a 
second-order slip boundary condition. Further, the high-Knudsen asymptotic expression 
for the flow rate is: 

Qoo ~ ^ log(Kn) + 0(l), (33) 

which implies a slow logarithmic divergence of the flow rate as a function of the Knudsen 
number. These two asymptotic limits ensure that the flow rate must have at least one 
minimum at some flnite Kn. The low-Knudsen asymptotic behavior is related to the 
situation where on the average the number of collisions encountered by a molecule in the 
bulk is much larger than the collisions with the wall. Thus, the effective balance between 
the frictional forces (due to colhsions) and the applied pressure gradient ensures a parabolic 
velocity proflle (with a slip at the wall). On the other hand, any molecular motion in high 
Knudsen number flows is retarded mostly by collisions occurring at the wall. This leads to 
an effectively flat velocity proflle in the channel. 

In addition to the asymptotic analysis, more detailed investigations of the linearized 
Boltzmann equation are available (see, for example, [2H])- In the next subsection, we 
compare the result of our numerical simulation with the numerical result of j2Hl and the 
asymptotic expressions. 
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5.1 Numerical simulation 

Poiseuille flow in a two-dimensional duct was simulated with the entropic LBGK model 
for a range of Knudsen numbers. The length to height ratio used in the simulation was 30, 
and the resolution was taken 1101 x 45 points. In lattice units, the Knudsen number is: 



The simulation results for the flow-rate is plotted in Fig. |21 together with the numerical 
solution of the stationary linearized BGK model |2H] and the asymptotic expressions. Some 
of the important conclusions following from this study are: 

• Quantitative agreement with the fully microscopic results is obtained in the range of 
< Kn < 0.01. This indicates that the parameter-free ELBGK model can be used 
in the domain of slip-flow for quantitative simulations. 

• At higher Knudsen numbers, the simulation results predict the expected logarithmic 
divergence (as indicated by the dashed line in Fig. refFig:Flowrate), but they do 
deviate from the numerical solution of the linearized stationary BGK equation. At 
very high Knudsen numbers, the finite size of the duct will always lead to a nearly 
flat (independent of Kn) fiow rate. 

• Our simulations predict the Knudsen minimum at moderate values of Kn. The mini- 
mum is found around Kn = 0.35, whereas much more elaborate (and computationally 
expensive) molecular dynamics simulations j2j| and the semi-analytical solutions to 
the stationary BGK equation indicate the location of the minimum at Kn in the 
range 0.8 — 1. The discrepancy is not surprising because the model used in our sim- 
ulation is drastically simpler than any of the alternative models. In particular, the 
built-in isothermal assumption limits its domain of validity to the cases when incom- 
pressibility is a good approximation indeed (slip- fiow). In the domain of Kn ~ 1 it 
is therefore required to use the non-isothermal model which includes correctly the 
energy redistribution in the collision processes. 

6 Discussion and conclusions 

In this paper, we have set up the basis of a new computational approach for the simulation 
of microfiows - the entropic lattice Boltzmann method. We have described the hierarchy 
of entropic lattice Boltzmann models for both isothermal and thermal simulations. Here, 
we have focused on the simplest isothermal model - the entropic lattice BGK equation 
and presented in detail the implementation of the diffusive wall boundary condition in 
the lattice Boltzmann setting. For the ELBGK model, we derived the parameter map of 
two-dimensional isothermal simulations on the three-dimensional data, and have tested all 
this on the benchmark problem of the planar Poiseuille flow in the full range of rarefaction 




(34) 
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Figure 2: Comparison of the ELBM solution for the dimensionless flow rate with the 
low Kn asymptotic solution, and the numerical solution of the stationary linearized BGK 
equation [2HI- The dashed line indicates the qualitative agreement with the expected 
logarithmic scaling at higher Kn. All the curves become indistinguishable at Kn < 0.01. 

(Knudsen number) from the nearly-continuous case to the free-molecular flow. Thus, we 
confirm the validity of the entropic lattice Boltzmann method as a viable tool for com- 
putations of the microfiows in the most relevant to MEMS applications domain of Kn up 
to Kn ~ 0.1. This confirmation clearly points to the usefulness of development of more 
sophisticated lattice Boltzmann models in this range of parameters which is the subject 
of further studies. It should be stressed that applications of the isothermal models can be 
considered as a parameter-free slip-flow hydrodynamic models, and that the way to extend 
the domain of the quantitative predictions must take into account energy conservation in 
the collisions. We would also like to point out the computational efficiency of the proposed 
models. For example, the full data set presented in Fig. El was computed within several 
hours on a single-processor computer facility. 

Finally, let us briefly mention the following by-product of our study, relevant for the 
use of the lattice Boltzmann simulations in molecular models. In particular, the Direct 
Simulation Monte Carlo method (DSMC) requires a good initial choice of the velocity 
distribution function to perform efficiently. Most of the current simulations use the local 
equilibrium (or equilibrium) Maxwell distribution function. A better choice could be the 
anisotropic Gaussian approximation 

f{v, X) ~ exp [-{Va - Ua)pP~p{Vp - Ufs)) , (35) 
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where P ^ is the inverse of the pressure tensor, 



Paf3 = / 7: fdv. (36) 



The LB data for the stationary stress tensor and for the flow field can be used as initial 
guess for the functions Ua{x) and P~p{x) in the velocity distribution function (|^. This 
initial condition for the DSMC simulation of the stationary will dramatically reduce the 
simulation time as compared to the case when it is initialized at equilibrium. 
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